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ABSTRACT 

We have succeeded in obtaining magnetized star models that have extremely strong 
magnetic fields in the interior of the stars. In our formulation, arbitrary functions of 
the magnetic flux function appear in the expression of the current density. By appro- 
priately choosing the functional form for one of the arbitrary functions which corre- 
sponds to the distribution of the toroidal current density, we have obtained configu- 
rations with magnetic field distributions that are highly localized within the central 
part and near the magnetic axis region. The absolute values of the central magnetic 
fields are stronger than those of the surface region by two orders of magnitude. By 
applying our results to magnetars, the internal magnetic poloidal fields could be 10 17 
G, although the surface magnetic fields are about 10 15 G in the case of magnetars. 
For white dwarfs, the internal magnetic poloidal fields could be 10 12 G, when the 
surface magnetic fields are 10 9 — 10 10 G . 
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1 INTRODUCTION 

The magnetic field inside a star is scarcely detectable by di- 
rect observations but has been considered to affect stellar evo- 
lutions and activities in many aspects. For instance, if strong 
magnetic fields are hidden inside degenerate stars such as 
white dwarfs or neutron stars, they may significantly affect the 
cooling process of the stars by providing an energy reservoir or 
by modifying heat conduction. Highly localized, anisotropic 
and relatively strong magnetic field configurations, on the 
other hand, may affect accretion modes onto degenerate stars 
in close binary systems by providing a well-focused channel 
of accretion to their magnetic poles. In order to know the pos- 
sible distributions and strengths of the magnetic fields inside 
the stars, we have to rely on theoretical studies. Until very re- 
cently, however, theoretical investigations could give us few 
hints about the interior magnetic fields. The reason for that 
may be twofold: one is related to the difficulty of the evolu- 
tionary computations of stellar magnetic fields and the other is 
related to the lack of methods to obtain stationary configura- 
tions of the magnetized stars. 

Concerning the evolution of the stellar magnetic fields, 
since it has been very difficult to pursue evolutionary com- 
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putations of the global magnetic fields for both interiors and 
exteriors of stars, few results have been obtained. Recently, 
however, Braithwaite and his collaborators have succeeded 



(Braithwaite & Spruit 


2004; 


Braithwaite & Nordlund 


2006; 


Braithwaite & Spruit 2006; Braithwaite 2006; Braithwaite 


2007; Braithwaite 200! 


ilBraithwaitdl2009tlDuez et alj|2010). 



They found that the twisted-torus configurations of the mag- 
netic fields inside stars seem to be stable across the dynamical 
timescale. 

On the other hand, to investigate possible structures of 
the interior and exterior magnetic fields by imposing station- 
arity is a different theoretical approach. Concerning this prob- 
lem, many attempts have been made but it has also been diffi- 
cult to obtain stellar structures with both poloidal and toroidal 
non force-free magnetic fields self-consistently, not only in 
the Newtonian gravity but als o in general re l ativity (see e.g. 
[ Chandrasekhar & Fermil 1953 ; Ferraro 19541 Chandrasekharl 
Il95d ' 



Prendergast 1956 



1960; Wentzel 1961 



Chandrasekhar & Prendergas^fl956 
Woltierlll959al; IWoltjerJll959bl: IWoltiei 
Ostriker & Hart wickll 19681; lMiketinac||l973l: iMiketinadl 19751 
Bocauet et al.l 1 19951: TToka & Sasakil |2004|; Kiuchi & Yoshidal 
Haskell et alj|2008l ; iDuez & Mafhisll2010h . It is only 



200: 



recently that axisymmetric and stationary barotropic stel- 
lar structures have been successfully solved for configu- 
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rations with both poloidal and toroidal magnetic compo- 
nents (Tomimura & Eriguchi 2005; Yoshida & Eriguchi 2006; 
lYoshida et alj|2006l ; lLander & Jonesll2009l : lOtani et al.N2009f) 
in a non-perturbative manner. 

It should be noted that the twisted-torus magnetic con- 
fig uration that appears during t he evolutionary computations 
by iBraithwaite & Spruitl d2004l) is qualitatively the same as 
one of th e exact axisymmetric and stationary solutions ob- 
tained in lYoshida et alj {2006). Moreover, stable configura- 
tions of stellar m agnetic fields must h ave a twisted-torus struc- 
ture according to Braithwaite (2009). Concerning the stability 
analysis, this type of configurations is expected to be stable, 
while magnetic fields with purely poloidal configurations or 
purely to roidal configurations have been shown to be unstable 
(see e.g. lTavledll973| IWrighdll973l: iMarkev & Tavlerlll973t 
iFlowers & Rudermanfl 19771) . 

In this paper, we apply the formulation develope d by 
Tomimura & Eriguchil(l2005h . I Yoshida & Eriguchi] (l2006h and 



Yoshida et al. ( 2006) in order to find out how strong and lo- 



calized poloidal magnetic fields can exist inside stars, as far 
as equilibrium configurations are concerned. In this formula- 
tion, the electric current density consists of several terms with 
different physical significances which contain arbitrary func- 
tionals of the magnetic flux function. These arbitrary func- 
tions correspond to the degrees of freedom in magnetized 
equilibria. One of the arbitrary functionals in the expression 
for the electric current density corresponds to the current in 
the toroidal direction. By choosing this functional form prop- 
erly, we would be able to obtain equilibrium configurations 
of axisymmetric barotropic stars with highly localized and ex- 
tremely strong poloidal magnetic fields. 



2 FORMULATION AND NUMERICAL METHOD 

Since we employ the f ormulation develop ed by 

iTomimura & Eriguchil d2005l) . lYoshida & Eriguch I d2006l) . 
and lYoshida et alj d2006h 7~here we summarize the main 
scheme briefly and explain newly introduced parts in detail. 



(x) No electric current is assumed in the vacuum region. 

(xi) The barotropic equation of state is assumed : 

p = p(p). (1) 

Here p and p are the pressure and the mass density, respec- 
tively. Assumptions of axisymmetric and equatorial symme- 
tries as well as rigid rotation are adopted here in order to sim- 
plify our investigations. 

In a rotating star under a radiative equilibrium, there ap- 
pears to be meridional flow in special cases. However, we ne- 
glect it because the time scale is many order s of magnitud e 
larger than the (magneto)hydrodynamic one dTassoufe oO). 
Also, there is a suggestion that gradual diffusi on of the in- 
ternal magnetic fields drives a meridional flow dUrpin & Ravi 
1994). The time scale of this, again, is much larger than the 
(magneto)hydrodynamic one. Thus it is also neglected. 

Under these assumptions, the basic equations are written 
as follows. The continuity equation is expressed as 

V ■ ( P v) = , (2) 

where v is the fluid velocity. The equations of motion in the 
stationary state are written as: 



-Vp = -V<t> B + m 2 e R + -( J -xH 
p p\c 



(3) 



where 4> g , Q,, j, c and H are gravitational potential, angu- 
lar velocity, electric current density, speed of light and mag- 
netic field, respectively. Here we use the cylindrical coordi- 
nates (R, ip, z) and e_R is the unit vector in the indirection. 
The gravitational potential satisfies Poisson equation: 



A(f> g = AixGp , 



(4) 



where G is the gravitational constant. Maxwell's equations are 
written as, 



V • E = 4ivp e 
V-JT = 0, 



V x E = , 



(5) 



(6) 
(7) 



2.1 Assumptions and basic equations 

We make the following assumptions for the magnetized stars. 

d 

(i) The system is in a stationary state, i.e. — = . 

at 

(ii) When stars are rotating and have magnetic fields, the 
rotational axis and the magnetic axis coincide. 

(iii) The rotation is rigid. 

(iv) The configurations are axisymmetric about the mag- 
netic or the rotational axis, i.e. — — = 0, where we use the 

dip 

spherical coordinates (r, 6, p). 

(v) The configurations are symmetric with respect to the 
equator. 

(vi) There are no meridional flows. 

(vii) The star is self-gravitating. 

(viii) The systems are treated in the framework of non- 
relativistic physics. 

(ix) The conductivity of the stellar matter is infinite, i.e. the 
ideal magnetohydrodynamics (MHD) approximation is em- 
ployed. 



V x H = 4tt^ , (8) 

c 

where p e and E are the electric charge density and the elec- 
tric field, respectively. Notice that we neglect the displacement 
current term in Eq.(|7} as is common in MHD approximation. 
The ideal MHD condition, or the generalized Ohm's equation, 
can be expressed as: 

E= -xH. (9) 

c 

We choose two kinds of barotropic equations of state. One is 
the polytropic equation of state: 



K p 1+1 ' N 



(10) 



where TV and Ko are the polytropic index and the polytropic 
constant, respectively. The other is the degenerated Fermi gas 
at zero temperature, defined as 



p = a[x(2x 2 - 3)\/x 2 + 1 + 31n(a;+ \J x 2 + 1)], (11) 
where 

P = bx, (12) 
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a = 6.00 x 10 dyn/cm 2 
b = 9.825 x 10Ve g/cm 3 



(13) 



(14) 



Here p e is the mean molecular weight. We fix = 2 in all 
our computations here, which corresponds to a fully ionized 
p ure hydrogen ga s. This choice of parameters is same as that 
m lHachisul Tl986). 



2.2 The form of the current density and the boundary 
condition 

From the assumptions of axisymmetry and stationarity, we in- 
troduce magnetic flux function ^ as follows: 



1 gjg 

7? dz ' 



H* = 



1 

RdR ' 



(15) 



where Hr and ff z are magnetic field components in the in- 
direction and 2-direction, respectively. We assume this flux 
function is positive in the entire space. By introducing this 
magnetic flux function, equation ^ can be automatically sat- 
isfied. It should be noted that the magnetic flux function can 
be expressed as: 



* = r sin 0A„ 



(16) 



where A v is the (^-component of the vector potential A — 
(A R ,A V ,A Z ). 

As shown in iTomimura & Eriguchi] d2005l) . for axisym- 
metric and stationary barotropes with rigid rotation we can 
constrain the form of the electric current density by using an 
integrability condition of the equations of motion, equation 
©= 

3 1 *«(¥) 



47T 



H + r sin 9 pLi($i)e lfi , 



(17) 



where «(\I r ) and fi(9) are arbitrary functions of the magnetic 
flux function Notice, in particular, that the toroidal compo- 
nent of magnetic field is given as 



— 



«(*) 



(18) 



which can be derived from equations ((3}, ([§) and ( 117) . It 
should be noted that these two arbitrary functions are con- 
served along the poloidal magnetic field lines. Although the 
meanings of the se two functions are described in previous 
works (see, e.g., lLovelaceetaL| [l986). in this paper we will 
explain their meanings differently from our point of view. 

Since we have assumed that there is no electric current 
in the vacuum region, in other words that there is no toroidal 
magnetic field outside the star (see equation [T7J, the form for 
k needs to be a special one. The simplest form can be k = 
constant dloka & Sasakil2004lHaskell et alj20o"3) . but for this 
choice of k the toroidal magnetic field would extend to the 
vacuum region. In order to avoid this possibility, we choose 
the functional form of n as follows: 



. 



(* - *max) 



fc + 1 



for 
for 



(19) 



This c hoic e of k is the same as t hat in lYoshida & Eriguchil 
J2OO6I) and lLander& Jonesl d2009h . In this paper we fix k = 
0.1. Equations d 1 8b and j 1 9b ensure that the toroidal magnetic 



field vanishes smoothly at the stellar surface. Incidentally, us- 
ing these functionals, we obtain the first integral of equation 
{3} as follows: 



dp 
P 



+ i(rsinf5) 2 fig + J M (*)d* + C\ 



(20) 



where C is an integration constant. The first term of the right- 
hand side is the gravitational potential. The second term on 
the right hand side is related to rotation. We can consider it 
as a rotational potential. Similarly, the third term means the 
potential of Lorentz force. We can regard this term as the 
magnetic force potential. Therefore, j /j, d^ is considered to 
be non-force-free contribution from the current density, as is 
seen in equation l !17| >. Since the Lorentz force is given by the 
cross product j/c x H, the first term of equation ( 117) has 
no effect on the equation of motion, i.e, it is force-free, and 
only the second term contributes to the Lorentz force, i.e., non 
force-free. The distribution of Lorentz force could be changed 
by adoptin g different functional form s for /x. All previ- 
ous works (Tomimura & Eriguchi 2005, Yoshida & Eriguchi 
| 2006|, lYoshida et alj|2006l lLander & Jonesl 120091 lOtani et al, 
120091) fixed /i = /io (constant). We choose a different func- 
tional form for [i in this paper as follows: 



Mo(* + c) m , 



771 + 1 



(* + <■) 



(21) 
(22) 



where m and e are two constant parameters. In order to avoid 
singular behavior, we fix e = 1.0 x 10 -6 in all calcula- 
tions. As we shall see below, the parameter m determines a 
degree of localization of the interior poloidal magnetic field. 
We assume that poloidal magnetic fields extend throughout 
the whole space and that there are no discontinuities even at 
the stellar surface. The global magnetic field configurations of 
our models are nearly dipole-like because of the requirement 
of the functional form for k at the stellar surface. These con- 
figurations contain closed poloidal magnetic field lines inside 
the star. The flux function \& attains its maximum at the cen- 
tral parts of these closed field lines and it takes its minimum 
on the symmetric axis and at infinity. The minimum value 
is zero because of ^ = r sin OA^ and the boundary condi- 
tion for A v — at infinity. The magnetic potential (J /id^ 
changes its qualitative behavior in its spatial distribution when 
m = —1. If we adopt m < — 1, as \& decreases from its max- 
imum to zero on the axis of the star the value of the magnetic 
potential increases unboundedly if e — > 0. As a result, the 
poloidal magnetic field lines are concentrated near the axis in 
order to fulfill such magnetic potential distributions. On the 
other hand, if we choose m > — 1, the value of the magnetic 
potential decreases as ^ decreases from its maximum to zero, 
which is realized on the axis. Then the poloidal magnetic field 
lines are distributed more uniformly than those for configura- 
tions with m < —1. If we choose m = 0, we obtain p = con- 
stant configurations. They are the same as those investigated 
by other authors. 

It is remarkable that the only freedom that we can take in 
our formulation is related to the choices of functional forms 
and the values of the parameters which appear in those func- 
tions. It implies that degrees of freedom for choices for these 
functions and parameters correspond to degrees of freedom 
for many kinds of stationary axisymmetric magnetic field con- 



© 0000 RAS, MNRAS 000, 000-000 



4 K. Fujisawa, S'i. Yoshida and Y. Eriguchi 



figurations. In fact, as we see from our results, different val- 
ues for m result in qualitatively different distributions for the 
magnetic potentials and the poloidal magnetic fields. In other 
words, we can control the magnetic field distributions to a cer- 
tain extent by adjusting the value for m. This is the reason why 
we use this functional form of /j, in this paper. 

After we choose the functional form of the current den- 
sity, by using Eq.J 17t and the definition of the vector potential, 
we obtain the following partial differential equation of the el- 
liptic type: 



for polytropic models and 

U = J g(x)d s r, 

g(x) = a{8x 3 [(x 2 + l)i - 1}} - p, 
for the Fermi gas configurations (see. Chandrasekhar 1939b , 



(33) 



(34) 



n = 



(35) 



A^sinv?) = A-kSa^, 8)sinp , Sa = - — • 



(23) 



As we have seen in the previous paragraph, all the physi- 
cal quantities related to the vector potential can be expressed 
solely by 4'. Therefore we need not solve for Ar and A z . It 
implies that our present formulation does not depend on the 
gauge condition for the vector potential A. Next we impose 
the boundary conditions for the gravitational potential and the 
vector potential, chosen as follows: 



A, 



O 



. r 



(r — > oo) 



This boundary condition for A v results in 



oi4 



(r 



(24) 



(25) 



(26) 



where H v is the poloidal magnetic field. From these boundary 
conditions and using a proper Green's function for the Lapla- 
cian, we have the integral representations of equation l(4j and 
equation J23b as follows: 



Mr) 



-G 



A v (r) s'mp — — 



Sa(v') sirup' d z r , 



(27) 



(28) 



Therefore, we can obtain smooth potentials, cj> g and A v by in- 
tegrating these equations. Since we have chosen the functional 
form of the current density which decreases near the surface 
and vanishes at the stellar surface sufficiently smoothly, we 
obtain continuous poloidal magnetic fields from A v , 



2.3 Global characteristics of equilibria 

To see the global characteristic of magnetized equilibria, we 
define some integrated quantities as follows: 



W=\J 4> aP d 3 r, 
T=\j p(m) 2 d 3 r , 



IT 



pd r , 



U = NIJ, 



(29) 



(30) 



(31) 



(32) 



(36) 



K — y (V x A)- Ad :i r = j H ■ Ad 3 r , 

where W, T, II, U, H and K are the gravitational energy, ro- 
tational energy, total pressure, internal energy, magnetic field 
energy and magnetic helicity, respectively. In order to evalu- 
ate the structures of magnetic fields, we define some physical 
quantities related to the magnetic fields as follows: 



H sur 



Clo r ' 2 M sine\H(r s ,e)\d6d<p 



s 



(37) 



where r s (8) and |-ff S ur| are the stellar radius in the direction 
of 9 and the surface magnetic field strength, respectively, and 
the surface area of the star is defined as: 



S = 



r s (8) sin 8 d8 dip. 



(38) 



The volume-averaged magnetic field strength in the central re- 
gion of the star is defined as 



C£ Jo" r 2 sin8\H(r,8)\drd8dp 



V 



(39) 



where we choose r c = 0.01r e and V is the volume of the 
central region with r < r c , defined as 

(•2tt 



V = 



r sin 8dr dddp . 



(40) 



This central region seems to be very small, but we can resolve 
it sufficiently because we use non-uniform and centrally con- 
centrated meshes (see. Fig. IA1 l and Eq. lA21l in Appendix). We 
have 77 meshes to resolve the region in actual numerical com- 
putations. 

In order to know the contributions of the poloidal mag- 
netic field and the toroidal magnetic field separately, we define 
the poloidal magnetic energy H p and the toroidal magnetic 
energy Ht as 



1 

8^ 



r 2 sin8\H r (r,8) 2 + Hg(r,8) 2 \drd8d<p , (41) 



Ht = 



-/ / / r 2 sm8\H v {r,8) 2 \drddd<p . 
if Jo io Jo 



(42) 



As for the magnetic multipole moment seen outside a star, we 
compute each multipole component by solving the following 
equation in a vacuum: 



A (A v sin ip) = . 



(43) 



Considering the boundary conditions at infinity and the sym- 
metry of the magnetized stars, the solution of the above equa- 
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tion can be expressed as 

oo oo 

A v wup = ^2 A,p, n sin tp = ^ b nA r~ n ~ 1 Y nil (6, ip) , (44) 

n=l n— 1 

where Y n ,i(9, tp) is the spherical harmonics of degree n and 
order m = 1. The coefficients b n .i correspond to the magnetic 
multipoles. 

2.4 Setting for Numerical Computations 

For numerical computations, the physical quantities are trans- 
formed into dimensionless ones using the maximum density 
Pmax, the maximum pressure p ma x and the equatorial radius 
r e as follows: 



r 



I 1 Pmaj 



for polytropic configurations and 



a 6 47rGp^ 



(45) 



(46) 



(47) 



for the Fermi gas models, and 

Pmax 

Here a is introduced so as to make the distance from the cen- 
ter to the equatorial surface of the star to be unity. Arbitrary 
functions are also transformed into dimensionless ones. Quan- 
tities witrTare dimensionless. For example, the dimensionless 
length is r and the dimensionless arbitrary functions are fi and 
k, respectively. Dimensionless forms of other quantities are 
collected in Appendix |ATJ 

The computational domain is defined as < 6 < ^ 
in the angular direction and < f < 2 in the radial direc- 
tion. Since the equation of magnetohydrostationary equilib- 
rium is defined only inside the star and the source terms of the 
elliptic equations for the gravitational potential and the mag- 
netic flux function vanish outside the star, our computational 
domain covers a region of the space that is sufficient for ob- 
taining equilibria. In order to resolve the region near the axis 
sufficiently, we use a special coordinate in actual numerical 
computations. Total mesh numbers in r-direction and in 0- 
direction are 1025 and 1025, respectively. We describe details 
of the computational grid points in Appendix IA2I 



2.5 Numerical method 



We use the scheme of iTomirnura & Erig uchi (2005.(). This 
scheme i s based on the Hachisu Self-Consistent Field (HSCF) 
scheme (Hachisu 1986), which is the method for obtaining 
equilibrium configurations of rotating stars. We define the ra- 
tio of the equatorial radius to the polar radius as the axis ra- 
tio q. This quantity q characterizes how distorted the stars are 
due to non-spherical forces. The stronger the non-spherical 
force becomes, the more distorted the stellar shape is. The 
non-spherical force can be the centrifugal force, the magnetic 
force or both of them. We fix the value of q in order to obtain 
the magnetized equilibria. We also fix one of /to and Clo- If we 
fix /to, we will obtain the value of Clo after the relaxation and 



VC 

-2.5 

-3 

-3.5 

_ "4 
O 

^ -4.5 
ai 
o 

-5 
-5.5 

-6 
-6.5 

1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 
log(Nr) 

Figure 1. The virial quantity VC, plotted against the number of grid 
points in the r-direction. 



iteration. If we fix f2o, we will obtain fio- Then, we will obtain 
one magnetized equilibrium state. 

2.6 Numerical accuracy check 

In order to check the accuracy of converged solutions, we com- 
pute a relative value of the virial relation as follows: 

I2T + w + m + u\ 



VC = 



\w\ 



(48) 



Since this quantity VC must vanish for exact equilibrium con- 
figurations, we can check the global accuraci es of the numer - 
ically obtained models as a whole (see e.g. lHachisulll986T) . 
Since the numerical results depend on mesh size, we have 
computed the same model by changing the number of grid 
points in the r-coordinate but fixing the number of grid points 
in the ^-direction as ng = 513. Fig.[TJshows VC as a function 
of the number of grid points in the r-coordinate for polytropic 
models. Since we use schemes of second-order accuracy, VC 
decrease s as the square inverse of the number o f grid points 
(see also Lander & Jones 2009; Otani et al. 2009). 



3 NUMERICAL RESULTS 

We give a brief summary of our numerical results here. First 
we show the basic features for negative m models and the de- 
pendences of the magnetic field configurations on the values 
of m for barotropes. We also show rotating and magnetized 
polytropic models in order to examine the effect of rotation 
on magnetic fields. The influence of the equation of state on 
the interior magnetic field is also displayed. We have com- 
puted N = 0.5, 1, 1.5 polytropic models and four white dwarf 
models with p c = 1.0 x 10 7 , 1.0 x 10 s , 1.0 x 10 9 , and 
gem 



3.1 Effect of the distribution of the toroidal current 
density on the distribution of the magnetic field 

We show the results for the distributions of the magnetic fields 
for different values of m. In particular, in order to examine 
the effect of magnetic fields alone, we consider configurations 
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Figure 2. Contours for the magnetic flux function (left panels) and for the logarithm of the strength of the magnetic field normalized 
by the averaged surface magnetic field (right panels) are shown. The inner solid circle corresponds to the surface of the star and the 
outer solid circle denotes the boundary of our computational region. The flux difference between two adjacent contours of thick lines 
is 1/10 of the maximum value of In the left panels, the thick poloidal field line is the boundary of the toroidal magnetic field region. 
The toroidal magnetic field exists only inside the region. In the right panels, the distribution of the logarithm of the magnetic field 
normalized by the averaged surface magnetic field, log 1( ) | H / H sur | contour, is shown. The thick solid curve corresponds to the curve 
with log 10 \H/H sur \ = 0. Inside this curve log 10 \H/H aur \ > and outside this curve log 10 \H/H sur \ < 0. The difference 
between two adjacent contours is 0.2. 
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5.25E-3 


8.31E-7 


2.789E-6 


0.5 


7.0E-5 


5.69E+0 


9.963E-1 


1.83E-5 


3.33E-1 


5.07E-2 


1.92E-1 


8.23E-7 


2.789E-6 


1.0 


6.1E-5 


4.78E+0 


9.959E-1 


1.70E-5 


3.33E-1 


5.07E-2 


6.90E+0 


8.11E-7 


2.789E-6 



Table 1. Physical quantities for firj = 0, kq = 10 and H Bur = 1.5 X 10 15 G polytropes with different values of m. 
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Figure 3. Profiles of log f3 at 9 = n/2, where ft is the plasma /3. Solid line represents the distribution for an m = 1.0 configuration, 
dashed line that for an m = 0.0 configuration and dotted line that for an m = —2.0 model, respectively. 



without rotation. The effect of stellar rotation is discussed in 
Sec. 13.21 Thus we set Qo — and compute TV = 1 polytropic 
equilibrium models with different values of m and appropriate 
values of q so that the surface magnetic field becomes roughly 
H sur = 10 15 G when p c = 1.0 x 10 15 gcm~ 3 and mass 
M = 1.4Mq. By setting N = 1 and an appropriate choice 
of polytropic constant K of p = Kp 2 , we obtain models with 
M = 1.4Mq. It should be noted that these models have the 
typical mass and radius for neutron stars. We choose TV = 1 as 
a simple approximation of neutron stars here. We searched and 
found the value of q by calculating many equilibrium states. 

Physical quantities of these models are shown in TableQ] 
It can be seen that values of II/| W| and a are almost the same 
among these models. Although the strength of the averaged 
surface magnetic field is H sur = 1.5 x 10 15 G, the values 
of "H/\W\ are much smaller than those of II/ \W\. It implies 
that the effect of the magnetic fields in these configurations 
on their global structures is very small. On the other hand, 
values of H c /H sur and H/\W\ vary rather considerably for 
different values of m. As the value of m is decreased, values 
of H c /H 3ur and H/\W\ increase. In Fig.|2]the structure and 
strength of magnetic fields are shown for three different values 
of to, i.e. to = —2.0 (negative to model), to = 0.0 (fi = con- 
stant model) and m — 1.0 (positive m model). The left-hand 



panels show the poloidal magnetic field lines and the regions 
where the toroidal magnetic field exists. The right-habd panels 
display the strength of the magnetic field |JT| normalized by 
the averaged surface magnetic field H aur . 

As seen from these figures, there are no discontinuities 
of the magnetic fields at the stellar surfaces. Due to the choice 
of the functional form of the arbitrary function and the 

distribution of the magnetic flux function, toroidal magnetic 
fields appear only in the region that is bounded by the outer- 
most closed poloidal magnetic field line inside the star (thick 
line). Thus the toroidal magnetic fields exist inside the torus 
region. 

As the value of to is increased, i.e. from top panels to bot- 
tom panels, the ratio of H c /H aur decreases (see left panels in 
Fig. O because the poloidal magnetic field becomes weaker. 
This is also related to the fact that the interior poloidal mag- 
netic field lines are much more localized near the axis for neg- 
ative to models. The contours of magnetic field strength also 
display the same tendency. For the to = 1.0 model, the con- 
tour of\H\ = H aU r (thick line) shows the stellar surface and 
the shapes of contours are nearly spherical. By contrast, the 
contours of the to = —2.0 model are highly distorted near 
the axis. The strength of the poloidal magnetic fields for the 
negative m models could exceed 10 17 G near the central re- 
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m = -1.5 m = 0.0 m = 1.0 




Figure 4. Isocontours of j v for different values of m. These panels show N = 1, q = 0.99 polytropic equilibrium models. The 
outermost curve denotes the stellar surface. The difference between two adjacent contours is 1/10 times the maximum of j v . The 
current is non-zero in the whole star except at the stellar surface and the symmetric axis (z-axis). 




gion. Fig.[3]shows the profiles of the plasma ft on the 9 — tv/2 
plane, i.e. on the equatorial plane. Here, the plasma j3 is de- 
fined as follows: 

P = 8tv P /\H\ 2 . (49) 

This quantity denotes the contribution of the gas pressure ef- 
fect compared with the magnetic pressure effect. Fig.|3]shows 
profiles of p for models with m = —2.0,0.0, 1.0. As seen 
from Fig|3] the profiles of /3 are very similar to each other 
near the stellar surface regions. For the region around f ~ 0.6, 
however, the value of /3 for the m = —2.0 model is larger than 
those for the m — and m = 1.0 models. Since these models 
have almost the same mass density distributions, this differ- 
ence means a difference of magnetic pressure distribution. In 
this region the magnetic field of the m = —2.0 configuration 
is weaker and thus the /3 becomes larger. However, it should 
be noted that these contours for the model with m — —2.0 are 
rather confined to the very narrow region near the central part. 
In other words, the gradient of the magnetic field distribution 
for the model with m = —2 is much steeper than the gradient 
of the gas pressure distribution compared with the models with 
m = 1.0 and m = 0.0. Thus the value of log ft becomes dra- 
matically small within the f [0 : 0.1] region and the minimum 
value of P can reach about ~ 20 in the central part. Therefore, 
in the central region of the model with m = —2.0 the influ- 
ence of magnetic field on the local structure of the star is no 
longer negligible. 

Here we explain the reason why this kind of highly lo- 



calized poloidal magnetic field configuration can be realized. 
We need to note the distribution of the toroidal current density 
j v in order to analyse our models properly, because the cur- 
rent density is related to the magnetic field closely by the two 
equations (8) and In Fig. [4] we show the distributions of 
the toroidal current density for models with different values of 
m. As seen from Fig. [4] the distribution of the toroidal current 
density is concentrated toward the magnetic axis for the con- 
figuration with negative values of m. This is due to the depen- 
dence of /} on the value of m. The current density distribution 
spreads over a large region inside the star as the value of m 
increases (from left panel to right panel). In other words, the 
distribution of the magnetic flux function becomes more and 
more concentrated toward the magnetic axis as the value of m 
decreases. It implies that the strengths of magnetic fields for 
models with negative values of m become very great near the 
magnetic axis. Our results show one possibility that a strong 
poloidal magnetic field can exist deep inside a star. If such a 
strong poloidal magnetic field is sustained deep inside a star, 
the contours of the magnetic field strength are no longer nearly 
spherical as in the bottom right panel of Fig. [2] Although this 
feature might be modified by dropping the assumption of the 
axisymmetry, it would give us one possibility for the presence 
of a strong poloidal magnetic field configuration deep inside a 
star. 

Finally, to characterize the magnetic structure we show 
the magnetic multipole moments of magnetized stars. In Fig. 
|5]the values of \b n ,i/bi,i | (equationl44t are plotted. There ap- 



© 0000 RAS, MNRAS 000, 000-000 



Axisymmetric and stationary magnetized barotropic stars 9 



pear to be only multipolar magnetic moments with odd de- 
gree (n = 1, 3, 5), because we have assumed the equatorial 
symmetry. As seen from these figures, in configurations with 
negative values of m the higher order magnetic multipole mo- 
ments contribute (|i>n,i/bi,i|) significantly to the total mag- 
netic field, while in configurations with positive values of m 
the magnetic dipole moment is the dominant component of the 
total magnetic field. These figures show that the external mag- 
netic field is nearly dipole when we adopt m — but it is not 
simple dipole when m > and m < 0. From the left panel, 
we see that the n — 3 (octupole) component reaches about a 
few tens of per cent of the dipole component when m = —2.0. 



3.2 Effect of stellar rotation 

We calculate two sequences with rotation for different values 
of m in order to examine the influence of rotation. We choose 
the value of p,o by obtaining a configuration with f2o = 
and q — 0.99 as a non-rotating limit of our equilibrium se- 
quence. We choose q = 0.99 here for simplicity. The value of 
q — 0.99 corresponds to an equilibrium configuration with 
H sur ~ 10 15 G when we consider a typical neutron star 
model with negative m. We have obtained sequences of sta- 
tionary configurations by fixing the parameters m and /to and 
changing the value of q. By changing the value of q for a fixed 
value of fro, we have equilibrium configurations with shapes 
that are deformed from spheres by rotational effect in addi- 
tion to the magnetic force. Since we fix the magnetic poten- 
tial parameter /to and m along one sequence, the equilibrium 
sequence is the one with approximately constant magnetic ef- 
fect. If the values of m and /to are changed, we will be able 
to solve another stationary sequence. We have calculated two 
stationary sequences with negative m (m = —1.5) and with 
m = 0.0, i.e. p, — constant . 

Physical quantities of stationary configurations are tabu- 
lated in Table [2] As seen from this table, the quantities jVl^j 
and a or the ratio n/| W| and T/| W| depend on the strength 
of the rotation. By contrast, magnetic quantities are almost un- 
affected by rotation. The dependence of the ratio H c /H sur 
on rotation is relatively small. The equilibrium configurations 
with highly localized magnetic fields that we have obtained in 
this paper are almost unchanged even by rapid rotation. There- 
fore, we do not consider the effect of rotation any longer in this 
paper. 



3.3 Effect of equations of state 

Thus far, we have discussed our magnetized configurations 
by showing the results for N = 1 polytropic models. The 
distribution of the toroidal current density, however, depends 
on the mass density profile through equation ( 117) . Thus we 
show other polytropic models, i.e. TV = 0.5 and N = 1.5 
polytropes, as well as configurations for degenerate gases, i.e. 
white dwarf models, in order to examine the influence of equa- 
tions of state on configurations with highly localized magnetic 
fields. 

We set q = 0.99 for polytropes and q = 0.999 for 
degenerate gases. The degenerate model with q — 0.999 
corresponds to a configuration with a H sur ~ 1.0 x 10 9 G 



magnetized white dwarf with m = —3.0, the central den- 
sity is 1.0 x 10 8 gcm -3 . This central density results in 
a white dwarf of about I.I6M0. Neither models rotates. 
We calculate 1 1 models with fixed values for q by setting m = 
-3.0, -2.5, -2.0, -1.5, -1.1, -0.9, -0.5, 0.0, 0.5, 1.0, 1, 3 
and examine the dependence of H c /H aur on the equation of 
state. 

Fig.|6]displays the ratio H c /H sur against the value of m 
for different equations of state. The dependency of this ratio 
on the value of m is qualitatively similar for these equations 
of state. Whichever equation of state we choose, we obtain 
configurations with highly localized magnetic fields, for which 
H c /H sur can exceed 100. The same is true for white dwarfs 
with highly localized magnetic fields. However, H c /H sur 
tends to become smaller for stiffer equations of state, as seen 
from Fig. [6] 

Fig. |7] and Fig. [8] display the distribution of mass den- 
sity, current density and the contour of log 10 \H\/H sur of 
m = —0.99 configurations. Fig. [7] shows results for poly- 
tropes N = 0.5 and N = 1.5 (stiffest and softest equations of 
state among the polytropic models considered here) and Fig. 
[8]shows results for white dwarfs with p c — 1.0 x 10 7 gcm~ 3 
and p c — 1.0 x 10 10 gcm~ 3 (stiffest and softest among the 
white dwarf models considered here). As seen from top pan- 
els in each figure, the mass density distributions of the softer 
equation of state (N = 1.5 and p c = 1.0 x 10 10 gcm -3 ) are 
more centrally concentrated than those of the stiffer equation 
of state (N = 0.5 and p c = 1.0 x 10 7 gcm~ 3 ). The cur- 
rent density distributions are also more centrally concentrated 
compared with the mass density distribution (middle panels). 
As a result, the poloidal magnetic fields become more highly 
localized for the softer equation of state (bottom panels). The 
mass of the white dwarf becomes higher for the higher central 
density. This implies that higher mass white dwarfs can have 
stronger interior magnetic fields deep inside if the magnetic 
field structure is fixed as in the present study. 



4 DISCUSSION AND CONCLUSIONS 

In this paper we have constructed axisymmetric and station- 
ary magnetized barotropes that have extremely strong poloidal 
magnetic fields around the central region near the magnetic 
axis. The strength of the magnetic field in that region could 
be two orders of magnitude larger than that of the surface 
magnetic field. In the context of the neutron star physics, this 
would imply that there might be magnetars whose interior 
magnetic fields amounting to 10 17 G if we assume the surface 
field to be order of 10 15 G and that there might be magnetized 
white dwarfs with interior magnetic fields that reach 10 12 G 
when the mass is nearly the Chandrasekhar limit and the sur- 
face field is of the order of 10 9 G. 

Moreover, it should be noted that highly localized mag- 
netized stars could have higher order magnetic multipole mo- 
ments in addition to the dipole moment. Although in most 
astrophysical situations magnetic dipole fields have been as- 
sumed, we may need to consider configurations with contribu- 
tions from higher multipole magnetic moments for some sit- 
uations. In those cases, configurations with negative values of 
m might be used to analyze such systems. 
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\W\ 


H/\W\ 


n/\w\ 


T/\W\ 


a 


6 2 


k 


VC 










m = — 1.5 


Ao — 


5.070E-7 










0.99 


4.75E+1 


0.9979 


9.71E-2 


1.14E-4 


3.33E-1 


0.00E+0 


5.06E-2 


0.00E+0 


3.33E-6 


1.80E-6 


0.9 


4.54E+1 


0.9972 


8.0OE-2 


1.26E-4 


3.19E-1 


2.14E-2 


4.51E-2 


1.17E-2 


3.61E-6 


3.41E-5 


0.8 


4.36E+1 


0.9960 


6.16E-2 


1.45E-4 


3.02E-1 


4.67E-2 


3.87E-2 


2.36E-2 


3.95E-6 


1.34E-5 


0.7 


4.19E+1 


0.9940 


4.40E-2 


1.74E-4 


2.85E-1 


7.30E-2 


3.22E-2 


3.37E-2 


4.32E-6 


5.70E-6 










m = 0.0 


Ao = 


5.520E-2 










0.99 


6.98E+0 


0.9947 


9.60E-2 


2.31E-3 


3.33E-1 


0.00E+0 


5.03E-2 


O.OOE+O 


1.20E-4 


1.85E-6 


0.9 


6.63E+0 


0.9934 


7.90E-2 


2.35E-3 


3.18E-1 


2.14E-2 


4.48E-2 


1.16E-2 


1.15E-4 


1.24E-6 


0.8 


6.29E+0 


0.9913 


6.05E-2 


2.39E-3 


3.01E-1 


4.73E-2 


3.83E-2 


2.37E-2 


1.06E-4 


1.36E-6 


0.7 


5.99E+0 


0.9878 


4.32E-2 


2.39E-3 


2.83E-1 


7.40E-2 


3.18E-2 


3.37E-2 


9.24E-5 


1.51E-6 



Table 2. Physical quantities of two sequences with m = 0.0 and m = —1.5. 
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Figure 6. The value of H c / -Hsur is plotted against the value of m for polytropes (q = 0.99) and white dwarfs (q = 0.999). 



4.1 Higher order magnetic multipole moments with even 

n 

It should be noted that in the analysis of this paper only higher 
magnetic multipole moments with odd n — 21 + 1 where i 
is an integer, i.e. 2 2e+1 moments, appear and that there are no 
higher magnetic multipole moments with even n = 21. This is 
due to the choice of the current density. Our choice of the arbi- 
trary function and the assumption of the symmetry of * 
about the equator necessarily result in magnetic field distribu- 
tions that are symmetric about the equator. It implies that the 
magnetic field should penetrate the equator and that 2 2e type 
distributions that are confined the upper or lower half of the 
space of the equator are excluded. In order to obtain closed 
magnetic field distributions in the half plane above or below 
the equator, the current density must be chosen so as to flow 
in opposite directions above and below the equatorial plane. 
It also implies that we need to set the current density on the 
equator in the (^-direction to vanish. 

Concerning 2 2i mu ltipole magnetic moments distribu- 
tions, ICiolfi et alj J2009I) have obtained such configurations. 
Their solutions correspond to the choice of the current density 
distributions that are antisymmetric about the equator. 



4.2 Forms of arbitrary functions 

One might think it curious that functions appear in the for- 
mulation and that there is no physical principle specifying 
how to choose those arbitrary functional forms. The same 
situation appears for the problem of calculating equilibrium 
structures or stationary structures of rotating and axisym met- 
ric barotropes. For that problem, the three component equa- 
tions of the equations of motion do not remain independent but 
come to depend on each other. This implies that one could not 
solve for all the three components of the flow velocity com- 
pletely. Assumptions of the stationarity and barotropy reduce 
the problem to a degenerate problem concerning the compo- 
nents of the flow velocity. Although there are three component 
equations for the three components of the flow velocity, those 
three component equations are no more independent. They be- 
come dependent each other due to the barotropic nature of the 
assumption for the gas. Therefore, one needs to specify the ro- 
tation law or corresponding relation in order to find stationary 
or equlibrium configurations for axisymmetric barotropes. The 
form of the rotation law is arbitrary. 

The only requirement for the functional form regarding 
the rotation law comes from the nature of the stability of the 
system. However, one needs to know the stability of the system 
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N = 0.5 N=1.5 




0.5 1 1.5 2 0.5 1 1.5 2 

R R 

Figure 7. Isocontours for p and ^ (top), j v (middle) and log[|f/|/|H sur |] (bottom). These panels are for configurations with 
m = —0.99. The left panels are contours for the model of a N = 0.5 polytrope and the right panels are for a N = 1.5 polytrope. 
The difference between two adjacent contours is 1/10 times the maximum of the corresponding quantities. 
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Figure 8. Same as Fig[7]except for the equations of state. The left panels are for a white dwarf with 1.0 X 10 7 gcm A and the right 
panels are for a white dwarf with 1.0 X 10 10 gcm -3 . 
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beforehand. If one does not have any information about the 
system to be solved, one has no principle by which to choose 
the form of the rotation law. 

The situation is the same for the stationary problem for 
axisymmetric magnetized bawtropes. For the stationary states 
of axisymmetric magnetized barotropes, the situation is more 
complicated than that for rotating barotropes, because not only 
the flow velocity but also the magnetic field appears in the 
problem. That also leads to the appearance of a greater num- 
ber of arbitrary functions in the problem. Thus it is very hard 
to specify the forms of arbitrary functions physically mean- 
ingfully. In such situations the only thing one can might be to 
explore many kinds of arbitrary functions to find out the gen- 
eral consequences of the resulting magnetic fields. 

Of course, if one could obtain a lot of information of the 
magnetic characteristics about the equilibrium states at hand, 
one could constrain the arbitrary functions more appropriately 
and more physically meaningfully. One possibility is to rely 
on the stability nature of the equilibrium, as in the rotating 
barotropic stars. Since there is no useful stability criterion for 
the field configuration with both poloidal and toroidal fields 
and linear stability analysis of the equilibrium is beyond our 
scope, we leave this issue of constraining the functional form 
for a future study. 



4.3 Application to magnetars 

The typical strength of the surface magnetic field of anoma- 
lous X-ray pulsar (AXP) and soft gamma-ray repeater (SGR) 
is considered to be 10 14 — 10 15 G by assuming the mag- 
netic dipole spin down (see e.g. iKouveliotou et all 1 19981; 
Kouveliotou et al. 1999; Murakami et al. 1999; Esposito et al. 
20091 ; lEnoto et alJl2009l ; TEnoto et al.ll2010l) . According to re- 
cent observational evidences, some types of AXP and SGR 
are regarded as similar kinds of isolated neutron star and 
are categorized as magnetars, although they were first con- 
side red to belong to two different types of neutron star, (see 
e.g. iDuncan & Thompson! 1 19921; iDuncan & Thompson! 1 199(j ; 
IWoods & Thompsonll200fj ; lMereghettill2008l) . 

For neutron stars with a strong magnetic field, such 
as magnetars, the strength of the maximum toroidal mag- 
neti c field inside has been e s timated to be 10 17 G (see 
e.g.lThompson & Duncanll 19951 ; iKluzniak & RudermaiJI 19981 ; 
Spruit 1999; Spruit 2009). Many authors have considered that 
only toroidal magnetic fields could become extremely strong 
and be hidden below the surfaces of the stars. Concerning 
poloidal magnetic fields, a very strong field is not consid- 
ered because it would be observed as a strong surface field 
since it is dipole-dominated. However, as shown in this pa- 
per, extremely strong poloidal magnetic fields can exist in 
the very central region at r c ~ 0.01r e , as seen from Tables 
[Tj and [2] and Fig. [6] and the definition of H c , Equationl|39t. 
If we apply our equilibrium models with negative values of 
m to magnetars with mass 1.4Mq, central density p max = 
1.0 x 10 15 gcm and average strength of the surface mag- 
netic fields 10 15 G, the strengths of the poloidal magnetic 
fields could be 10 16 — 10 17 G. We also consider weak mag- 
netized magnetars w ith average stren gth of the surface mag- 
netic fields 10 13 G teea et al]|2010h . If we apply our equi- 
librium models, the strengths of the poloidal magnetic fields 



could be 10 14 - 10 15 G. Since these strong poloidal magnetic 
fields located nearly along the magnetic axis in the central core 
region, the magnetic structures in the core region are highly 
anisotropic. If extremely strong magnetic poloidal fields are 
hidden within the core region, there could be magnetic fields 
with higher order multipole moments. 

If the neutron star shape is deformed by a strong mag- 
netic field and the magnetic axis is not align ed the rota- 
tional axis, gravitational waves will be e mitted dCutleril2002l ; 
lHaskell et al]|2008l ; lMastrano et aT]|201 lh . Gravitational wave 
emission tends to become stronger as the ellipticity of the 
meridional plane of the star becomes larger. For our models, 
decreasing m increases the value of 1 — q in the H sur constant 
sequence (see the value of 1 — q in TableQJ. Thus those mod- 
els with highly localized magnetic field here may be efficient 
emitters of gravitational wave. 



4.4 Some features of highly magnetized white dwarfs 

It is widely believed that the effect of the stellar magnetic 
fields play a significant role in astrophysics. For example, iso- 
lated magnetized white dwarfs tend to have a higher mass 
than n on-magnetic white dwarfs fwickramasinghe & Ferrariol 
2000). According to observations, the surface magnetic field 
strength of white dwarfs varies f rom very little to 10 9 G 
( Wickramasinghe & Ferrario 2000). Therefore, there are some 
strongly magnetized white dwarfs whose surface magn etic 
field about 10 8 -10 9 G. For example, Ijordan et al.l l l 1998b es- 
timated the field range 3.0 x 10 8 -7.0 x 10 s G in GD 299. 
EUVE J03 17-855 is a massive high-field magnetic white 
dwarf with rapid rotation. Its magnetic field was calculated 
by an offset dipole model with 4.5 x 10 8 G and period of 
725 s. PG 1031+234 is a high-field magneti zed white dwarf. 
ISchmidt et all dl986h and lLatteretakl l l 19871) estimated its ro- 
tation period 3.4 h and its magnetic field as 5.0 x 10 8 - 
1.0 x 10 9 G. The observed spectral variations cannot be fit- 
ted well by a simple dipole magnetic or offset dipole model, 
so they have proposed a two-component model composed of 
a nearly centered dipole and a strongly off-centered dipole. In 
other words, the magnetic field structures of several strongly 
magnetized white dwarfs could not be explained by applying 
simple dipole structures. 

We have obtained strongly magnetized white dwarfs with 
higher order magnetic multipole moments in this paper. If we 
apply our configurations with negative m, some strongly mag- 
netized star such as PG 1031+234 may have strong interior 
magnetic fields. According to our numerical results, H c could 
reach as high as 10 12 G when H 3ur ~ 3.0 x 10 9 G for a 
highly localized (rn = —3.0) and high mass (p c — 1.0 x 10 9 , 
M ~ 1.34M ) model (see Fig. HJ. Since the central mag- 
netic field strength H c depends on the equation of state as we 
have shown in Sec l3.3l it becomes higher as the central den- 
sity increases. Thus high mass white dwarfs could have strong 
poloidal magnetic fields according to our models with negative 
m. As we have displayed in Sec. 13. II N = 1.5 polytropes with 
negative values of m have rather large higher order magnetic 
multipole moments. The same is the case for magnetized white 
dwarf models, i.e. they have rather large higher order magnetic 
multipole moments. Therefore, the magnetic fields outside of 
such stars are far from simple dipole fields if the magnetized 
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white dwarfs have highly localized strong poloidal magnetic 
fields deep inside the stars. 



4.5 Comments on stability of magnetized barotropes 

Once equilibrium configurations are obtained, it would be de- 
sirable to investigate their stability. However, a satisfactory 
formulation for the linear stability analysis for general mag- 
netic configurations has not been fully developed, although 
there is a stabilit y criterion o nly for purely toroidal magnetic 
configurations ( ITavledll973l) . For purely poloidal or mixed 
poloidal-toroidal magnetic configurations, magnetic configu- 
rations with rotation or other general situations, no authors 
hav e ever succeeded in obtaining a clear sta b ility criterion (see 
e.g. IMarkev & Tavlej|l973l; Iwrighdll973t. iMarkev & TavleJ 
1 1974 lTavleJll98dl ; iBonanno & Urpinll2008|) . Therefore the 
stability of the configurations obtained in this paper contain 
both poloidal and toroidal magnetic fields has not been inves- 
tigated. 

By contrast, the stability of magnetized stars may be 
investigated through that time-dependent evolutionary com- 
putations of the magnetic configurations. Thanks to power- 
ful computers, some authors have recently employed mag- 
netohydrodynamical codes to follow the time evolutions of 
magnetized configurations and find out whether these con- 
figurations would settle down to certain 'stable equilibrium 
states'. Such investigations concerning the magnetic config- 
urations have been carried out by Braith waite and his cowork- 
ers as mentioned in Introduction (see e.g. Braithwaite & Spruit 



2004 ; IBraithwaite & Nordlundl 120061; iBraithwaite & Spruij 
20061; lBraithwaitejl2006l ; lBraithwaitell2007l ; lBraithwaitell2009l " 



Duez et al.ll2010h . According to their results, purely toroidal 



configurations and purely poloidal configurations are shown 
to be a ll unstable, as prev i ously s hown or expected (e.g. Tavler 
1 19731 ; IMarkev & TavTeil Il973l; Flowers & Rude rmanl Il977 



However, see 



Geppert & Rheinhardtj |2006 for some results 



about stability). Concerning the mixed poloidal-t oroidal mag- 
netic configurations, r ecent numerical studies dBraithwaitd 
l2009l ; lDuez et alj2010h have shown that they are stable as long 
as the following condition is satisfied: 



(50) 



where Qo is a numerical factor of 10 — 10 depending on the 
stellar structures. By performing 3D MHD simulations, it has 
been shown that non-axisymmetric perturbations to equilib- 
rium stars grow when this condition is not satisfied. Stars with 
mixed magnetic fields whose dominant component is poloidal 
field seem to evolve toward non-axisymmetric configurations 
until the amplitude of the perturbations reach nonlinear regime 
and saturate. As can be seen from tables in this paper, we have- 
found no models that satisfy that criterion (equation I50t for 
our particular choice of functional forms presented above (see 
Sec |2.2| >, because the energy stored in the toroidal magnetic 
field is at most a few per cent for all of our models. In order 
to obtain configurations that satisfy the criterion, we need to 
choose different functional forms from those used in this pa- 
per. We should be careful to apply the criterion, however, to 
general configurations of magnetic fields. The class of solu- 
tions with both toroidal and poloidal magnetic fields obtained 



here may be rather different from the ones studied by Braith- 
waite and his collaborators, even if they share the obvious 
characteristics of twisted-torus structures of magnetic fields. 
As is seen in completely diffe rent stability natures of seem- 
ingl y similar configura tions in iGeppert & Rheinhardtl d2006l) 
and IBraithwaite! ( |2007|) . it is quite uncertain at this moment 
that failure to satisfy the criterion (equation 150) for our mod- 
els here means unstable nature of them. It would be interesting 
to study the stability nature of our configurations thorough ei- 
ther linear perturbation analysis or direct MHD simulations. 

4.6 Conclusions 

In this paper, we have presented an extended formulation for 
obtaining axisymmetric and stationary barotropic configura- 
tions with both the poloidal and toroidal magnetic fields. We 
have shown the possibility that magnetized stars have strong 
poloidal magnetic fields inside the star. Our findings and con- 
jectures can be summarized as follows. 

(i) By choosing the functional form for one of the arbitrary 
functions that appear in the basic formulation for the configu- 
rations under the assumptions mentioned before, we have ob- 
tained magnetized configurations in which extremely strong 
poloidal fields are confined within the central part of the near 
axis region. When we apply our models to magnetars, the inte- 
rior magnetic strength would be 10 17 G while the surface mag- 
netic strength is 10 14 - 10 15 G. On the other hand, if we apply 
our models to magnetized white dwarfs with mass ~ 1.34M© , 
the surface field strength would be 10 9 G and H c reaches 10 12 
G. 

(ii) If stars have extremely strong poloidal magnetic fields 
deep inside, the contours of magnetic field strengths are not 
spherical but rather column-like shapes as shown in the fig- 
ures. 

(iii) If stars have extremely strong magnetic fields deep in- 
side, contributions from higher order magnetic multipole mo- 
ments to the outer fields around the stars cannot be neglected. 
This implies that if stars have highly localized and extremely 
strong magnetic fields deep inside, then observations of mag- 
netic fields around the stars could not be explained by the sim- 
ple dipole models that have been used in most situations. 
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APPENDIX A: NUMERICAL METHOD 
Al Dimensionlcss quantities 

In this paper, physical quantities are used in their dimension- 
less forms as follows: 



9 47rGr|p max 



Q. EE 



n 



K 



H = 



H. 



suffix 



H suf fix 



max 



K - 



v47rGr|p max 
K 

47rGrfpi iax ' 



(Al) 
(A2) 
(A3) 
(A4) 

(A5) 
(A6) 
(A7) 
(A8) 
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3v 



C 



C 



(A9) 



(A10) 



where m and 712 are the mesh numbers defined as follows: 



47rGrlp ma x 

Here H su ffix is the component of the magnetic field where 
suffix may be c (center), sur (surface), p {poloidal) and t 
(toroidal). Similarly we define normalized global quantities as 
follows: 

M 



M 



fepmax 



W = 



W 



T 



n 



u = 



u = 



47rGr|p2 

n 

47rGrfp2 
U 



H 



max 



(All) 



(A12) 



(A 13) 



(A14) 



(A 15) 



(A 16) 



We also define dimensionless forms of arbitrary functions as 
follows: 



0, 

«(*) = < ko 



k- 



\fc+l 



for * < 
for * > * n 



(A17) 



m = -(n r — 1) + 1 
n 2 = -(n r - 1) + 1 



(A22) 
(A23) 



Here n r is the total mesh number in the r-direction. In prac- 
tice, since we use a difference scheme of the second-order ac- 
curacy for the derivative and Simpson's integration formula, 
we divide each mesh interval defined above further into two 
equal size intervals in the r coordinate. We use n r — 513 and 
thus the actual total number of the mesh points is (2n r — 1) = 
1025. 

Concerning the ^-direction, we have to resolve the region 
near the axis, because for m < values the magnetic fields 
seem to be highly localized to the axis region. In order to treat 
such magnetic fields near the axis region, we introduce the 
following mesh in the ^-direction: 



= x j , x j = (J 



I) AX, l<j<n e , AA = 



ne — 1 

where ne is the total mesh number in the ^-direction. We also 
divide each mesh interval defined above further into two equal 
size intervals. Then, we use ne = 513 and thus the actual total 
number of the mesh points is 1025. Fig. lAll shows the relations 
between the order of the grid points and the r- or ^-coordinate 
value. 



(A24) 



dk(t>) _ J , for * < * n 

~ 1 ko(* - * mal )' , for * > * n 

for k and 



Mo 
m + 1 



(* + e) 



(A18) 

(A19) 
(A20) 



for jit. We choose k = 0.1, fto = 10 and e = 1.0 X 10 6 and 
keep their values fixed during all calculations in this paper. 



A2 Computational grids 

We describe the details of our numerical grid points. In or- 
der to resolve the distributions of the source term of the vec- 
tor potential equation without loss of accuracy, we choose 
the following non-uniformly distributed grid points in the ac- 
tual numerical computations. In the f-direction, we divide the 
whole space into two distinct regions: [0, 1.0] (region 1), and 
[1.0, 2.0] (region 2). In each region, the following mesh points 
are defined: 

Wi = (i — l)Aiui , Aiui = — ^5 , for 1 < i < m , 

' = w ? s J2-V1 ' (A21) 

uii — 1.0 + (i— ni)Au>2, Au>2 = , for m < i . 

«2 — 1 
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Figure Al. Left: the coordinate r/r e is plotted as a function of the grid points. The solid curve shows the region 1 ([0, 1]) and the 
dashed curve shows the region 2 ([1, 2]). Right: the same as the left panel except for the 9 coordinate. 
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